Modelo de regresión

\[y(x)=\eta(x;\theta)+\varepsilon\quad, \]

con \(\theta\) los parámetros desconocidos del modelo, de interés para el experimentador.

Diseño Experimental

• Diseño exacto: una colección de puntos, con sus réplicas, en los que tomar medidas.
• Diseño aproximado: un diseño, \(\xi\), visto como una medida de probabilidad.

FIM

La Matriz de Información de Fisher se define como

\[ M(\xi,\theta)=\sum_{x\in\mathcal{X}}f(x)f^t(x)\xi(x), \]

con \(f(x)=\partial\eta(x, \theta)/\partial\theta\)

Criterios de optimalidad

• D-optimalidad: \(\Phi_D=|M^{-1/k}|\).
• \(D_s\)-optimalidad: \(\Phi_{D_S}=\left \{\frac{|M|}{|M_{22}|}\right\}^{-1/s}\).
• A-optimalidad: \(\Phi_A=\text{Tr} (M^{-1})\).
• I-optimalidad: \(\Phi_I=\text{Tr} (AM^{-1})\), con \(A=\int_{\mathcal{R}}\mu(x)f(x)f(x)^Tdx\).

Teorema General de Equivalencia

Un diseño \(\xi^\star\) es \(\phi\)-óptimo si y sólo si,

\[Tr \left\{\nabla \phi(\xi^\star)[M(\xi^\star, \theta)-f(x)f^t(x)]\right\}\leq 0, x\in\mathcal{X}.\]

La igualdad se cumple en los puntos de soporte del diseño.

Eficiencia

La eficiencia de un diseño \(\xi\), para un cierto cirterio \(\Phi\), se define como

\[ eff_\Phi(\xi) = \frac{\Phi(\xi)}{\Phi(\xi^\star)} \]

Algoritmo Cocktail

1. Empezar con un diseño no singular.
2. Añadir \(x\) tal que \(\max_{x\in\chi} f^t(x)\nabla \phi(\xi)f(x)\).
3. Optimizar los pesos a través del algoritmo multiplicativo.
4. Aplicar heurísticos.
5. Comprobar optimalidad, si no, repetir 2-5.

Cálculo de diseños óptimos: función opt_des

resArr.D <- opt_des("D-Optimality",
  model = y ~ a*exp(-b/x),
  parameters = c("a", "b"),
  par_values = c(1, 1500),
  design_space = c(212, 422))

Output

kable(resArr.D$optdes)

resArr.D$sens

resArr.D$convergence

Distribuciones de probabilidad

Se puede seleccionar con el parámetro distribution

• ‘Homoscedasticity’.
• ‘Gamma’, que puede ser usada también para una exponencial o una normal heterocedástica con error relativo constante.
• ‘Poisson’.
• ‘Logistic’.

también se puede dar la función directamente con el parámetro weight_fun.

Diseño para la distribución de Poisson.

poissionDopt <- opt_des("D-Optimality",
   y ~ t0 + t1 * x + t2 * x^2, 
   c("t0","t1","t2"),
   c(0.95, -1, -1), 
   c(0, 0.5799), 
   distribution = "Poisson")

Heurísticos del algoritmo

FP2_01 <- opt_des("I-Optimality",
  y ~ a0 + a1*log(x) + a2*log(x)^2,
  c("a0", "a1", "a2"),
  c(1, 1, 1),
  c(0.0001, 1),
  reg_int = c(0.0001, 1))

que tiene una eficiencia de \(81.83\%\).

Heurísticos del algoritmo

FP2_02 <- opt_des("I-Optimality",
  y ~ a0 + a1*log(x) + a2*log(x)^2,
  c("a0", "a1", "a2"),
  c(1, 1, 1),
  c(0.0001, 1),
  reg_int = c(0.0001, 1),
  join_thresh = 0.005)

con una eficiencia mayor a \(99.99\%\).

Eficiencia de un diseño

Un diseño cualquiera

design <- data.frame("Point" = c(220, 240, 400),
  "Weight" = c(1/3, 1/3, 1/3))

design_efficiency(resArr.D, design)

[1] 0.3063763

O los dos diseños calculados antes

design_efficiency(FP2_02, FP2_01$optdes)

[1] 0.8183762

optdes

print(resArr.D)

     Point    Weight
1 329.2966 0.5000068
2 422.0000 0.4999932

summary(resArr.D)

Model: 
y ~ a * exp(-b/x)
and weight function: 
NULL
Optimal design for  D-Optimality :
     Point    Weight
1 329.2966 0.5000068
2 422.0000 0.4999932

 Minimum efficiency (Atwood):  99.9986396401789%
 Criterion value:  9972806

Clase optdes

plot(resArr.D)

Plot del diseño óptimo

Región de puntos candidatos

get_augment_region("D-Optimality",
               init_design = resArr.D$optdes, 
               alpha = 0.3,
               model = y ~ a * exp(-b/x),
               parameters = c("a", "b"),
               par_values = c(1, 1500),
               design_space = c(212, 422),
               calc_optimal_design = F)

Región de puntos candidatos

Diseños aumentados

aug_arr <- augment_design(
               "D-Optimality",
               init_design = resArr.D$optdes, 
               alpha = 0.3,
               model = y ~ a * exp(-b/x),
               parameters = c("a", "b"),
               par_values = c(1, 1500),
               design_space = c(212, 422),
               calc_optimal_design = F)

Eficiencia de \(80\%\)

Añadiendo los puntos 260 y 380 con peso 0.15

aug_arr

     Point    Weight
1 329.2966 0.3500048
2 422.0000 0.3499952
3 260.0000 0.1500000
4 380.0000 0.1500000

Algoritmo de redondeo

1. Empezar con un diseños aproximado con \(k\) puntos y pesos \(\xi(x_i)\)
2. Calcular \(\nu = n - k/2\)
3. Asignar \(z_i \gets \nu\xi(x_i)\)
4. Redondear \(z_i\) al entero más cercano \(r_i = \lceil {z} \rceil\)
5. Si \(\sum_i r_i = n\) entonces devolver \(r = \{r_1, \dots, r_n\}\) si no incrementar \(r_i\) con \(i = argmax_i (n \xi(x_i) - r_i)\) o disminuir \(r_i\) con \(i = argmin_i (n \xi(x_i) - r_i)\) hasta \(\sum_i r_i = n\)

Redondeo

design_test <- data.frame("Point" = seq(1, 5, length.out = 7),
  "Weight" = c(0.1, 0.0001, 0.2, 0.134, 0.073, 0.2111, 0.2818))

efficient_round(design_test, 21)

Comprobar eficiencia

exact_design <- efficient_round(resArr.D$optdes, 7)
aprox_design <- exact_design
aprox_design$Weight <-
  aprox_design$Weight / sum(aprox_design$Weight)

design_efficiency(resArr.D, aprox_design)

[1] 0.9897433

Shiny Antoine

shiny_optimal()

Shiny Addpoints

shiny_augment()

Resumen

• Facilidad de uso
• Aumentar diseños para diferentes modelos
• Redondeo eficiente para poder implementar los diseños

Futuras mejoras

• Modelos con varias variables explicativas
• Interfaz gráfica con funcionalidad completa

El paquete en github y CRAN

• https://github.com/Kezrael/optedr
• https://cran.r-project.org/web/packages/optedr/index.html

¡Gracias por vuestra atención!